Setup
Load the rwrfhydro package.
library("rwrfhydro")Set the path to the directory of WRF Hydro test cases.
fcPath <- '~/wrfHydroTestCases/Fourmile_Creek_testcase_v2.0'Visualize Gridded Domain Data: VisualizeDomain()
When working with a new model domain, it can be highly informative to explore the physical or geospatial nature of the model through visualization. The WRF-Hydro model has gridded geospatial information for both the LSM and hydrologic components, which are potentially at different spatial resolutions. The hydrologic components may also employ a separate file containing vector information, which we do not discuss here.
This section demonstrates a basic functionality for visualizing the gridded geospatial information in R. These functions suitable for medium sized domains. Proper GIS software is better suited for visualizing large domains and more intensive spatial representations of model components.
The same function VisualizeDomain helps visualize both LSM and hydro components. It plots the information over map data which can be retrieved from a variety of sources, with the help of the ggmap package.
To properly display the spatial information, R currently requires the use of the rgdal and sp packages though these are not required for general installation and use of rwrfhydro. We hope to relax this assumption in the near future. An underlying function GetDomainCoordsProj does the work of obtaining and transforming the projection and/or datum of the spatial data to render it a lat-lon WGS84 map.
Land surface model pixels
We begin by taking a look at the spatial information file for the LSM. First we specify the path to the ‘Fourmile_Creek’ test We pick the “geo” file, which holds the geospatial information for the LSM. The following commands plot the topographic height for each cell over a map.
geoFile <- paste0(fcPath,'/run.ChannelRouting/DOMAIN/geo_em_d01.Fourmile1km.nlcd11.nc')
coordsProj <- GetDomainCoordsProj(geoFile)
ClosureGeo <- VisualizeDomain(geoFile, plotVar='HGT_M', plot=FALSE, plotDf=coordsProj)
mapMargin <- .05*c(-1,1)
closureRtnGeo <-
ClosureGeo(zoom=11, pointsize=9,
grad=topo.colors(15), alpha=.6, maptype='terrain',
subsetRange=range(plotDf$value),
xlim=range(plotDf$long)+mapMargin*1.5,
ylim = range(plotDf$lat)+mapMargin)## Warning: bounding box given to google - spatial extent only approximate.## converting bounding box to center/zoom specification. (experimental)## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040369,-105.442619&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN
The step of getting coordsProj could be done internally to VisualizeDomain but in some cases (e.g. plotting many variables for a big domain) can speed up the visualization.
For illustration purposes, we set plot=FALSE in the call to VisualizeDomain to skip the default plot output. Typically this first, default plot serves as a basis for tweaking the options to the returned function, named ClosureGeo in this example, which help achieve the desired graphic representation.
Note carefully that VisualizeDomain returns a function, such a returned function is termed a closure: a function with data inside it. (Because our style convention is first-letter capitals for functions, ClosureGeo begins with a capital letter.) The arguments to this function allow most of its main attributes to be tweaked when called. Because the graphics in this function are handled using the ggplot2 package, which has a signifcant learning curve, we design these closures to give the user the 90% solution without th need to learn ggplot2. The return values of the closures give the data and the ggplot objects which can befurther manipulated as necessary.
Critically, zoom and pointsize depend on our domain and even the size of the graphics device. To work more directly with the data and the plot object, see the closureRtnGeo object. The returnComponents argument to the closure function gives even more fine grained output. We’ll illustrate working with the return object below.
Routing pixels
A similar set of function calls allows visualization of the hydrologic grids. We look at the gridded channel network file. This time we skip the call to GetDomainCoordsProj.
hydroFile <- paste0(fcPath,'/run.ChannelRouting/DOMAIN/Fulldom_hires_hydrofile.Fourmile100m.nc')
ClosureHydro <- VisualizeDomain(hydroFile, plotVar='CHANNELGRID', plot=FALSE)
closureRtnHydro <-
ClosureHydro(zoom=11, pointsize=1,
## can reference the internal plotDf (or other variables internal to the closure)
location=c(lon=mean(plotDf$long), lat=mean(plotDf$lat)), alpha=.2,
grad=c('white','blue'), maptype='terrain')## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040415,-105.442563&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN
We highlight use of several of the available arguments to the closure. To see the full set of arguments to the closure, simply invoke the print method on the closure: call the function without paraentheses. This will return the entire code and the arguments are shown at the top with their default values (e.g. argument=default). Reading the coode will reveal two variables which actually come from the VisualizeDomain function: bbox and plotDf. (Because the returned function encloses this data, it is refered to as a closure). Several of the arguments to the closure can reference these variables: location, xlim, ylim, subsetRange.
While the above command shows all the hydrologic grid pixels, many of them are missing values (-9999) because they are not on the gridded channel network. We use the subsetRange keyword to limit the range of values represented in the plot. We again use the xlim and ylim keywords to zoom in. We also limit the color scale to blue.
closureRtnHydro <-
ClosureHydro(zoom=14, pointsize=1,
## can reference the internal plotDf (or other variables internal to the closure)
location=c(lon=mean(plotDf$long), lat=mean(plotDf$lat)),
subsetRange=c(0),
grad='blue', maptype='terrain',
xlim=mean(plotDf$long)+mapMargin/2,
ylim =mean(plotDf$lat)+mapMargin/2/1.5)## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040415,-105.442563&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN## Warning: Removed 258 rows containing missing values (geom_point).
To illustrate use of the return values of the above closures, we show how to combine the two plots.
names(closureRtnGeo)## [1] "plotDf" "ggObj"names(closureRtnHydro)## [1] "plotDf" "ggObj"closureRtnGeo$ggObj +
geom_point(data=closureRtnHydro$plotDf,
aes(x=long, y=lat), size=.5, shape=21, fill='white', color='darkblue') +
ggtitle('Fourmile Creek, CO - Elevation and Stream Channel')
The returnComponents argument to the closures will return more fine-grained objects from the plot in caseu those are useful.
